import matplotlib.pyplot as plt
import numpy as np
def f(x):
    return np.exp(x[0]+3*x[1]-0.1)+np.exp(x[0]-3*x[1]-0.1)+np.exp(-x[0]-0.1)

def gradient(x):
    return np.array([np.exp(x[0] + 3 * x[1] - 0.1) + np.exp(x[0] - 3 * x[1] - 0.1) - np.exp(-x[0] - 0.1),
                     3*np.exp(x[0] + 3 * x[1] - 0.1) - 3*np.exp(x[0] - 3 * x[1] - 0.1)],dtype=np.float32)

def hessian(x):
    return np.array([[np.exp(x[0] + 3 * x[1] - 0.1)+ np.exp(x[0] - 3 * x[1] - 0.1) + np.exp(-x[0] - 0.1),
                      3 * np.exp(x[0] + 3 * x[1] - 0.1) - 3 * np.exp(x[0] - 3 * x[1] - 0.1)],
                     [3 * np.exp(x[0] + 3 * x[1] - 0.1) - 3 * np.exp(x[0] - 3 * x[1] - 0.1),
                      9 * np.exp(x[0] + 3 * x[1] - 0.1) + 9 * np.exp(x[0] - 3 * x[1] - 0.1)]],dtype=np.float32)

def direction(x):
    return - np.linalg.inv(hessian(x)) @ gradient(x)


alphas = np.linspace(0.1,0.45,10)
betas = np.linspace(0.1,0.9,10)
epsilon = 1e-7
count = np.zeros((len(alphas),len(betas)))
for i in range(len(alphas)):
    alpha = alphas[i]
    for j in range(len(betas)):
        cnt = 0
        beta = betas[j]
        x0 = [1,1]
        while True:
            d = direction(x0)
            lamda = gradient(x0) @ d
            if -0.5 * lamda < epsilon:
                break
            cnt += 1
            t = 1
            x1 = x0 + t * d
            while f(x1) > f(x0) + alpha * t * lamda:
                t = beta * t
                x1 = x0 + t * d
            x0 = x1
        count[i][j] = cnt
# 三维图
x,y=np.meshgrid(alphas,betas)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, count)
ax.set_xlabel('alpha')
ax.set_ylabel('beta')
ax.set_zlabel('count')
plt.show()